#!/usr/bin/python3
# by Sun Smallwhite <niasw@pku.edu.cn>(https://github.com/niasw) & Shi Lei <shijim@gmail.com>

import numpy
import matplotlib
matplotlib.use('Agg'); # plot and save, never showing

import matplotlib.pyplot as plt

def calcXvector(clussprs):
  '''
# for a paper p:
#  x(p)=1/sqrt(|#c(p)|)
# input:
#  clussprs: cluster seperators: [0,p_1,p_2,...,p_n,|#P|]
  '''
  xvec=numpy.zeros(clussprs[-1])+1.;
  for c in range(1,len(clussprs)):
   xvec[range(clussprs[c-1],clussprs[c])]=1./numpy.sqrt(clussprs[c]-clussprs[c-1]);
  return xvec;

def monitorXandErr(eigvec,clussprs):
  '''
# monitoring x vector and error function
  '''
  xvec=calcXvector(clussprs);
  reerr=1-xvec.dot(eigvec)**2/xvec.dot(xvec)/eigvec.dot(eigvec);
  return xvec,reerr;

def visualizeDivide(eigvec,clussprs):
  '''
# visualize dividing progress
  '''
  mappedxvec=calcXvector(clussprs);
  mappedeigv=eigvec;
  maponexvec=mappedxvec/numpy.sqrt(len(clussprs)-1);
  plt.clf();
  plt.plot(range(0,len(eigvec)),mappedeigv,color="green",linewidth=2,linestyle='--',label='eigenvector');
  plt.plot(range(0,len(eigvec)),maponexvec,color="blue", linewidth=2,linestyle='-',label='x vector');
  plt.grid();
  plt.xlabel('paper index (sorted by eigvec)');
  plt.ylabel('vector component');
  plt.legend(loc='upper left');
  plt.savefig('../out/xveig'+str(len(clussprs)-1)+'.png');
  # plt.show(); # don't forget to commentize "matplotlib.use('Agg');" in the beginning

def clussprs2clusters(clussprs,sortmap):
  '''
# generate clusters set from cluster seperator data
  '''
  def _idxfunc(x):
    return x[1];
  clusters=[list(set([_idxfunc(sortmap[p]) for p in range(clussprs[c],clussprs[c+1])])) for c in range(0,len(clussprs)-1)];
  return clusters;

def calcClusters(eigvec, clusnum=20, monitor=False):
  '''
# set paper clusters to make x approching eigv
# by Sun Smallwhite <niasw@pku.edu.cn>(https://github.com/niasw)

# input:
#  eigenvector
#  clusnum: recording start point, max number of clusters
#  monitor: whether visualize the progress
# output:
#  snapshot (last history state detail): clusters set [[paper indices],...]
  '''
  dimension=len(eigvec);
  curcnum=dimension;
  setcnum=clusnum;
  bestrcd=numpy.abs(clusnum-curcnum); # best record
  noimprovecount=0; # prevent loop
  clussprs=[];

  def _keyfunc(x):
   return x[0];
  sortmap=sorted([(eigvec[p],p) for p in range(0,len(eigvec))],key=_keyfunc); # (eigen_value,origin_order)
  sortbck=sorted([(sortmap[p][1],p) for p in range(0,len(eigvec))],key=_keyfunc); # (origin_order,eigen_order)
  mappedeigv=numpy.float64([e[0] for e in sortmap]);
  numpyeigvec=numpy.array(mappedeigv);
  eigsum=sum(numpyeigvec);

  if (monitor):
   print('1st try:'+'\n  now cluster number: '+str(dimension)+'\n  set cluster number: '+str(setcnum));
  if (dimension<clusnum):
   print('Setting Cluster Number is smaller than Current Cluster Number, reset to equal.');
   clusnum=dimension;
   setcnum=clusnum;
   curcnum=clusnum;
   clussprs=list(range(0,dimension+1));
  while (curcnum!=clusnum and noimprovecount<=7): # no improvement for 7 times => treat as loop
   setcnum=int(numpy.ceil(setcnum*numpy.power(clusnum/curcnum, 0.5))); # setting cluster number more smoothly
   #setcnum=int(numpy.ceil(setcnum*numpy.sqrt(clusnum/curcnum))); # setting cluster number
   curcnum=0; # current cluster number

   clussprs=[0,dimension]; # seperators
   objecteigv=[sortmap[p][0]*numpy.sqrt(setcnum) for p in range(0,dimension)]; # objective eigenvector
   cursor=dimension;

   while (len(clussprs)<=setcnum):
    # print('cluster number: '+str(len(clussprs)-1)+'/'+str(len(eigvec)));
    cluspnum=int(numpy.ceil(1/objecteigv[cursor-1]**2)); # number of papers a cluster contains
    if (cursor<=cluspnum):
	     cursor=0;
	     # print('WARNING: seperation reaches end ... theoretic error may increase.');
	     break; 
    newcluspnum = cluspnum + 0.1;
    maxitcnt = 1;
    curitcnt = 0;
    while(newcluspnum > cluspnum and curitcnt <= maxitcnt):    	
	    cluspnum=int(numpy.floor(newcluspnum));
	    if(curitcnt>=maxitcnt):
	    	break;
	    curitcnt+=1;  	    
	    offset=int(numpy.ceil((cluspnum+1)/2.));
	    if (cursor<=offset):
	     cursor=0;
	     # print('WARNING: seperation reaches end ... theoretic error may increase.');
	     break;
	    newcluspnum=int(numpy.ceil(1/objecteigv[cursor-offset]**2)); # recalculate iteration once
	    if (cursor<=newcluspnum):
	     cursor=0;
	     # print('WARNING: seperation reaches end ... theoretic error may increase.');
	     break; 	      		       	        	    	    	    	
    #print('curitcnt:'+str(curitcnt));
    if (cursor <= 0):
    	break; 
    clussprs.insert(1,cursor-cluspnum);
    cursor=cursor-cluspnum;
    if monitor:
      xvec,reerr=monitorXandErr(mappedeigv,clussprs);
      visualizeDivide(mappedeigv,clussprs);
      print('unparallel ratio = '+str(reerr));
   curcnum=len(clussprs)-1;
   if (numpy.abs(clusnum-curcnum)<bestrcd):
    bestrcd=numpy.abs(clusnum-curcnum);
    noimprovecount=0;
   else:
    noimprovecount+=1;    
   numvec=numpy.array([(clussprs[it+1]-clussprs[it]) for it in range(0,curcnum)]);
   sqrtnumvec=numpy.sqrt(numvec);
   eigclussum=numpy.array([sum(numpyeigvec[clussprs[it]:clussprs[it+1]]) for it in range(0,curcnum)]);
   flowdotproduct=sum(eigclussum/sqrtnumvec);
   if monitor:
     print('improve to #cluster=' + str(curcnum));
     print('list of cluster size:'+str(numvec));
  clusters=clussprs2clusters(clussprs,sortmap);
  return clusters;

